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Abstract 



We present a method to compute the magnetic moment of a bulk, finite- 
size, three-dimensional, anisotropic superconductor. Our numerically im- 
plemented perturbative procedure is based on a solution of the nonlinear 
Maxwell-London electrodynamic equations, where we include the nonlinear 
relation between current and gauge invariant velocity. The method exploits 
the small ratio of the finite penetration depths to the sample size. We show 
how to treat the open boundary conditions over an infinite domain and the 
continuity requirement at the interface. We demonstrate how our method 
substantially reduces the computational work required, and discuss its imple- 
mentation to an oblate spheroid. The numerical solution is obtained from 
a finite-difference method. We briefiy discuss the relevance of this work to 
similar problems in other fields. 
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I. INTRODUCTION 



A large number of problems in electrodynamics and related areas such as fluid dynam- 
ics, involve the solution of partial differential equations for certain fields inside and outside 
a finite region of a definite geometrical shape. For many common geometries, and when 
the boundary conditions are simple (e.g. fields or their normal derivatives vanishing at the 
boundaries) the solution can be found, often with ease, from analytical or simple numerical 
methods. However, for more complicated situations where one has less trivial boundary 
conditions, or when the equations are made more complicated by the presence of nonlin- 
earities, analytical methods may be unavailable and numerical techniques encounter serious 
difficulties. 

One of these situations pertains to the electrodynamics of a superconducting sample 
of finite size. It is well-known that in the limit where the electromagnetic fields do not 
penetrate the sample the problem can be rather easily solved [ |[|. However, this is hardly 
ever the case of interest: the physical information one obtains in experiments comes in fact 
from the penetration of the fields inside the sample, characterized by penetration depths 
which, although small, cannot be neglected. 

Consider a superconductor that occupies a bounded, macroscopic region Q C R^, in the 
presence of an applied uniform magnetic field. Ha. For Ha below some critical value, super- 
conductors are in the so called Meissner regime [ ^, where the magnetic flux is expelled from 
the bulk of the sample. Their behavior is similar to that of material which is both an ideal 
conductor and an ideal diamagnet. The applied magnetic field generates a resistance-free 
current which produces a magnetic field that opposes Ha- As a consequence, everywhere 
except very close to the interface (within a few penetration depths), the magnetic field van- 
ishes: this is known as the Meissner effect. Except for the most trivial geometries such as 
infinite slabs or isotropic spheres, the relevant boundary value problem becomes then nu- 
merically very awkward: basically one is faced with solving the appropriate electrodynamic 
equations in the entire space, not just in Q, while the most important variation of the fields 
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takes place in a very thin boundary layer in Q. 

This question has come recently to the fore in the context of the study of high tem- 
perature superconductors (HTSC's). Identifying the symmetry of the paired electrons, the 
so called pairing state, in these materials, would provide an important clue to the still un- 
known mechanism responsible for superconductivity in HTSC's. It turns out [ that 
careful measurements of certain electrodynamic properties in superconductors, can provide 
fingerprints for the nodal structure of the order (gap) parameter [ |^ (i.e. the points or lines 
in the Fermi surface where it vanishes). The unconventional pairing states which are widely 
believed to exist in HTSC's, produce nonlinearities in the electromagnetic response. As we 
shall see, the resistance-free current, in addition to the usual terms linear in the superfluid 
velocity, includes a small contribution for which the current- velocity relation is nonlinear. 
These nonlinearities then give rise [ 0J^, in the Meissner regime, to a magnetic moment 
which has a small but detectable transverse component, m±, perpendicular to the field Ha 
even when this is applied along a direction of symmetry of the sample. This occurs when 
the applied field lies in the a — b crystallographic plane [ ^ (the z axis is taken to be along 
the c crystallographic direction). The angular dependence of the transverse component m±, 
as the crystal is rotated about the z axis, reflects directly the symmetry of the pairing state. 
It is this quantity that has been experimentally studied [ |] for purposes of the identification 
of the pairing state. 

The physics of the situation has been extensively discussed in [ |^ , and we deal here with 
the mathematical and numerical implications. The computation of the magnetic moment 
requires the solution of a problem of precisely the kind described in the previous paragraphs. 
One must solve the appropriate electrodynamics, the Maxwell-London equations described 
in Section |T|, for all space, since at infinity the boundary condition requires that H — > Ha. 
On Q these equations contain, as we shall see, important and nontrivial nonlinearities. A 
solution in the limit of zero penetration of the fields in the sample is possible, but it would 
be completely inadequate, since it would reflect only the geometry and not the detailed 
electromagnetic response of the superconductor. On the other hand, the numerical solution 
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for the nonlinear electrodynamics in all space would be computationally demanding. 

In this paper we present a discussion of the methods that we have developed to obtain 
results [ for the nonlinear magnetic moment, including both the longitudinal and trans- 
verse components, in HTSC's. These methods involve the numerical implementation of a 
perturbation scheme in which the small expansion parameter is the appropriately defined 
ratio of an effective penetration depth to a characteristic dimension of the superconductor. 
We will show that this numerical implementation reduces the problem essentially to that 
of finding the numerical solution in Q. In the region \ Q outside the sample, one turns 
out to need only a solution for the scalar Laplace equation with trivial Neumann bound- 
ary conditions. For a sufficiently symmetric Q, the form of the solution can be obtained 
analytically, while for some other numerical solution would suffice. 

In Section H we discuss the nonlinear Maxwell- London equations for a superconductor 
and show how they give rise to the magnetic moment. The geometrical shapes we have 
considered for Q (dictated by experimental considerations) are discussed in Section ^ where 
we introduce the general solution in R'^\i7. In Section |V| we discuss the computation of the 
magnetic moment and present the main result of this paper, the perturbative method and its 
numerical implementation. Numerical considerations for the equations in Q are described 
in the Section 0. In Section |V|, the equations are solved for the previously discussed 
geometry, using a modified Gauss-Seidel relaxation, with the nonlinear terms (which are 
nonanalytic) included through Picard's method. We also illustrate the general ideas of the 



perturbation method. Finally, in Section |VI1| we give conclusions and guidelines for possible 
improvements and generalizations. While our discussion is in the context of superconducting 
electrodynamics, it will not escape the reader's notice that our procedures can be adapted to 
many other problems in which one is faced with a small "skin depth" or similar parameter 
involving the penetration of fields or of their derivatives inside a region, and the problem 
can be more easily solved in the limit where this parameter vanishes. 
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II. MAXWELL-LONDON ELECTRODYNAMICS 



A. Maxwell-London Equations 

We begin by introducing the steady state Maxwell-London equations [ |^-|TT| that provide 
the framework to compute the field distributions we need in order to evaluate the magnetic 
moment. 

As stated in the Introduction we consider a superconductor in an applied uniform mag- 
netic field Ha that occupies a bounded simply connected region Q C and at its boundary, 
dQ, is surrounded by vacuum. On \ Q the current is j = and therefore in the steady 
state the local magnetic field H satisfies the Maxwell equations 

V-H = (2.1a) 

VxH = 0. (2.1b) 
The problem reduces to that of finding a magnetic scalar potential $ that satisfies 

H = -V$ (2.2a) 

= 0. (2.2b) 

On Q the relevant fields are H, the superconducting current j, and the "superfiuid velocity 
field" [ |l3 defined as: 



Vs = ^ + -A, (2.3) 
2 c 

where x is the phase of the superconducting order parameter, A the vector potential, and e 
the proton charge, (with h = = !)■ The field Vg is conventionally defined as above, and 
actually has units of momentum. The relation between Vg and H is given by the second 



London equation [[Ty]: 



V X Vg = -H. (2.4) 



In the steady state the appropriate Maxwell equation is Ampere's law, 



An 

V X H = — j, (2.5) 

c 



substituting Eq. ( |2.4| ) into (|2.5|) we obtain: 



V X V X = ^j(v,) (2.6) 



which is the general equation that will be investigated in this article. It must be supple- 
mented by a relation, j(vs), which will be discussed in |II C] , to be substituted in its right 



hand side (RHS). Eq. (p^ ) must then be solved together with Eq. ( p.2|) and the proper 
boundary conditions. 

The required boundary conditions are the following: first at infinity one must have, 

- V$ = Ha. (2.7) 



Second, deep inside the sample all fields must vanish. Third, H must be continuous [ |^,[Tl]|T2 
on the boundary dQ. Finally, the currents are confined to the superconducting material, 
j ■ n = in dQ, n is the unit normal pointing outwards. The first boundary condition (at 
an open boundary over an infinite domain) can be satisfied by construction of the solution 
in R'^ \ Q. The remaining boundary conditions have to be implemented in the numerical 
algorithm. 

Because of these boundary conditions we see that, as emphasized in the Introduction, 
this problem can indeed be computationally very demanding. It involves solving nonlinear 
differential equations, in principle in an unbounded region, but with the relevant fields 
varying very rapidly in a small region inside the material. We will see here and in the next 
Sections how these difficulties can be overcome, for suitable geometries, by making use of a 
numerical implementation of a perturbation scheme. 



B. The Magnetic Moment 



Let us at this point introduce some considerations that point to the eventual way out 
of these numerical difficulties. We recall that the quantity of interest here is the magnetic 

8 



moment m. From the current distribution j in Q, the magnetic moment can be obtained by 



volume integration [ |T3 



m 



Ij^rffir" xj(r"), (2.8) 



where r" is the position vector for a point in the region Q. 

An important consequence of nonhnear Maxwell- London electrodynamics and uncon- 
ventional pairing states is that m need not be aligned with the applied field Ha even if the 
latter is applied along a direction of geometrical symmetry, (along a principal axis of the 
demagnetization tensor [|l|] of the body). For simplicity, we will restrict ourselves to this case 
here, although it is straightforward to add the complications arising from a non-diagonal 
demagnetization tensor. 

Let us introduce here, for a typical macroscopic experimental HTSC sample, the ratio 
between some length d characterizing its size, and the characteristic value of the London 
penetration depth, which we shall denote by A. In a macroscopic sample the ratio e = X/d 
is a very small quantity. This small ratio is essentially what will be used in this work to 
develop a perturbation method to compute the magnetic moment to first order in e. The 
starting point of this procedure is the existence of a solution in the limit e ^ (when 
all the nonlinear effects vanish), which corresponds to imposing trivial Neumann boundary 
conditions on dQ. For a suitable choice of Q, such as an ellipsoid, the form of this solution 
may be found analytically. The perturbation method may also be applied, as we shall see, 
in certain cases when only a numerical solution in the small e limit is available. 

The components of m parallel and perpendicular to Ha apphed along a direction of 
symmetry, can generically be written for e <C 1 in the form 

m\\ = mo{l - a\\ e + 0{e^)) (2.9a) 

m^=mo(a±e+ 0{e^)), (2.9b) 

where tuq denotes the longitudinal magnetic moment in the limit A = 0, (and therefore 
e = 0), which is proportional to Ha- It depends only on the geometry of Q and therefore 



contains no physical information. For Q in the shape of an elhpsoid, values are given in [ 
|I[ . For finite A there is a reduction, linear in the field to leading order, in the absolute value 
of my. This reduction is due to current penetration in the material and it implies a positive 
constant a\\. For a very few simple geometrical shapes and linear, isotropic. Maxwell- London 
equations, values of a\\ are given in textbooks [ PJI^. When nonlinear effects are included, 
they contribute a correction to a\\ linear in the field, but their most conspicuous effect is the 
appearance, in general, of nonvanishing values of a±, also proportional to the field [ ||,||. It 
is for this reason that the transverse component is the physical quantity of interest. 

We next derive an alternative expression for m that helps to fully exploit the existence 
of the small parameter e. Using Eqs. ( |2.1| ), ( p.5|) and formulas for vector calculus [ ^ 
can transform the the quantity in the integrand of ( |2.8| ): 



one 



r X (V X H) = V(r-H) - (r- V)H-H (2.10) 



(r ■ V)H = -V X (r X H) - 2H. 



(2.11) 



After substitution of Eq. (^) into Eq. (|CT| ), 



r X (V X H) = V(r • H) + V X (r X H) + H, 



(2.12) 



integration over Q and use of Gauss' theorem yields: 

m = — / d^[n(r" ■H)+ nx (r" X H)] + — / c/^]H = ml + m2, (2.13) 
Stt Jan 8tc Jn 

where r" is the position vector for a point on dQ. The notation mi, m2 refers to the two 
terms in the middle portion of Eq. ( p.l3|) . If we recall Eq. ( |2.4| ) and use an alternative form 
of Gauss' theorem we can also rewrite m2 as a surface integral: 

mi = — / dS\n (r" ■ H) + n x (r" x H)], ma = — / dS n x v^. (2.14) 

Stt Jan Sne Jan 

The terms mi and m2 are of different order in e and the latter is small, i.e. of O(emo). 
This follows from the volume expression for m2, as seen in Eq. (|2.13| ): since H is confined 
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to a "skin" layer of thickness A from the surface, we see at once that m2 will be of order A 
times the applied field, that is, of order emg, thus explicitly vanishing in the zero penetration 
limit. Alternatively, from Eq. ( p.4| ) one sees that Vg scales as A and then the same result 
follows from Eq. ( p. 141 ). Specifically if the equations and boundary conditions require the 
field to decay exponentially far from the surface (up to polynomial corrections), then after 
decomposing the volume integral into surface and normal components we have 

X f r'^max _ JiL 3 

m2 = — dS' dw}_^{Hraax)ie < —Smax{Xi{Hmax)i)- (2.15) 

OTT JS' Jo ^ OTT 

This expression is not proportional to the volume of Q, V, as is the case for ttiq, but rather 
to \S ~ 0{eV) where S is the surface area of dQ. As a result m2 is O(emo). 



C. The relation between j and Vg 

We now return to the pending question of the equation relating the fields j and needed 
to supplement (|2.6|) . This is given by the usual two-fluid phenomenology [ ^,0: 



oo 



j(v,) = -eAr^/ d'sn{s)vf[{vf■^r,)+2 rfC /(i5^(C) + v/ ■ v,)] (2.16) 

FS J 

where Nf is the total density of states at the Fermi level, n{s) is the density of states at point 
s at the Fermi surface (FS), normalized to unity, v/(s) is the s-dependent Fermi velocity, / 
is the Fermi function, with E{() = ((^^ + |A(s)|^)^/^, T the absolute temperature and A(s) 
is the superconducting order (gap) parameter. 

In the simplest approximation (which we will call the linear case), one considers only 
the linear terms in the expression for j(vs), that is, one expands the right side of ( 2.16| ) and 
writes: 

where A is the penetration depth tensor. In the special case of an isotropic superconductor, 
and only in this case, the fields H, j and Vg all satisfy the vector Helmholtz equation. 
But HTSC's are in general highly anisotropic, layered structures with penetration depths 
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much smaller [ |^ in the layers {a — b planes) than in the direction perpendicular to them 
(along the c-axis), so that the isotropic limit does not apply. The components of A in its 
diagonal representation are the square of the London penetration depths, Xa, Xb and Ac and 
experimentally Ac ^ Xa, Xb- In the numerical examples discussed here we will, for simplicity, 
neglect the comparatively small in-plane anisotropy and use the notation Xab for the average 
of Aa and A;,. 

In the anisotropic case Eqs. (|2.9| ) are still valid with e expressed in terms of an effective 
penetration depth, primarily determined by whichever component Aj plays the dominant 
role in the current decay. For example, for the geometry considered in the next Section, the 
relevant quantity is the penetration depth in the crystallographic a — b plane. 

In the problem of interest here one must consider, instead of Eq. (|2.17|) , the nonlinear 
terms arising from the full relation (Eq. (|2.16|) ) between j, and Vg and substitute this in the 



RHS of Eq. ( p.6|) . After suitable assumptions for the FS and other physical quantities are 
introduced, this can be done by performing the FS integrals numerically, [ |^ or, in the low 
temperature limit, analytically [ Inclusion of these nonlinear terms is crucial, because, 
as mentioned, the physically important angular dependence of the transverse magnetic mo- 
ment arises precisely from these nonlinear effects. The actual expressions for j(vs) used here 
are taken from [ |^ and are quoted in Appendix A. 

III. ANALYTICAL CONSIDERATIONS: OBLATE SPHEROIDAL GEOMETRY 



Experimental samples [ |T6[ in which magnetic measurements are performed are in the 
shape of a flat "disk" with rounded edges, and the axis of revolution along the crystallo- 
graphic c axis of the crystal. This is well approximated by taking f2 to be a fiat ellipsoid 
of revolution (an oblate spheroid). For an oblate spheroid it is possible to find an analytic 
expression for the general form of the potential $. Therefore, our ideas can be implemented 
in this geometry in terms of analytic expressions in \ Q. 

The potential $, for e = 0, satisfies trivial Neumann boundary conditions on dQ, and the 
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solution contains a single parameter which is simply related to m||. When the penetration 
depth is finite, the longitudinal moment does change, but its correct value can in principle 
be determined from the boundary conditions and the solution inside, either through an 
iteration process as described in the next Section, or, much more efficiently, through the 
perturbation method we shall develop. 

The fundamental equation (|2.6|) is not separable in spheroidal coordinates. (Even the 
vector Helmholtz equation is not). Still, it is desirable to employ these coordinates as dQ 
is then described by a single parameter and this significantly simplifies the process of nu- 
merically fulfilling the boundary conditions. The simple implementation and discretization 
of the boundary conditions on dQ yields higher accuracy where it is most needed, since 
boundary grid points lie on the interface. 

We denote the major and minor semiaxes of the spheroid by A and C respectively, and 
we have A » C for actual samples. We take (see Fig. ||) a coordinate system fixed to 
the direction of the magnetic field, with its z-axis parallel to the c crystallographic direction 
of the superconductor (and parallel to the C semiaxis of an ellipsoid). The field is applied 
along the x-axis, and we picture the experiment as being performed by rotating the crystal 
about the 2;-axis and measuring the angular dependence of m_|_. As the crystal is rotated the 
axes X — y remain fixed in space, and should not be confused with the coordinates, affixed 
to the crystal structure, that we shall also use and denote by x', y'. We call the angle 
between and x'. 



In the definition we use [ |T^, the oblate spheroidal coordinates ^,ri,(p, are related to 



Cartesian coordinates by the transformation 

x = f{l + e')'/'(l - r/2)i/2 cos(^, (3.1a) 

y = f{l + eY^\l-v'Y^'smv, (3.1b) 

z = Kv, (3.1c) 
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where 0<^<oo,— l<?7<l,0<(y9< 27r, and / is a focal length scale factor (the 
coordinates ^, t] and ip are dimensionless) . In Fig. 2 we show this coordinate system 
at a fixed azimuthal angle if = 0°. One can obtain the relation between ^, rj, ip and 
Cartesian coordinates fixed to the crystal by replacing with + \1/ in For example 

x' = /(I + — r]^)^/^ cos((/9 + \E'). The relation between unit vectors in these and 

Cartesian coordinates is 

^ = (^2 /^2)i/2 (^(l - V'Y^' cosy^x + e(l - v'Y^' Sin y^y + (1 + e')'/'^^), (3.2a) 

^ = - (^2 /^2)i/2 ((l + ^'y^'v COS + (1 + ^2)1/2^ sin - ^(1 - r^')^/';^), (3.2b) 

= — sin + cos (3.2c) 

We see that ^ = n is the unit normal pointing outwards. In the limit / — > 0, the spheroidal 
system reduces to the spherical coordinate system. For / finite, the surface ^ = const 
becomes spherical as ^ — oo: 

f^^r, rj ^ cos 9, as ^ oo. (3.3) 

In the same limit i ^ r and r) — > ~9, where r and 9 are spherical coordinates. Various 
quantities in oblate spheroidal coordinates are given in Appendix B. 

To construct the general form of the solution for the fields in the region \Qwe employ 
an electrostatic analogy [[I|]. The current distribution is localized within Q and the magnetic 
potential in the exterior region can be written as an expansion in the appropriate set of 
orthogonal functions which in this case are the spheroidal harmonics [ (a generalization 

of spherical harmonics), characterized by angular and azimuthal indices I and m. Thus we 
write 

oo I 

^ = ^a + J2J2i^T cos(my,) + A- smim^))fQTmPriv) ^ > ^o, (3.4) 

1=0 m=0 

where Pf^ and Qf^ are the associated Legendre functions of the first and second kind respec- 
tively, and $a is the potential due to the applied field. The condition that — V$ HaX 
when ^ ^ oo yields 
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= -HafPlitOPiiv) cos<^. (3.5) 

Since for ^ ^ oo QY^{iC) ^0 \/ m, I this fulfills the boundary condition at infinity. The 
remaining terms in Eq. ( p.4|) are due to the presence of the superconductor. The coefficients 
and A\j_, for example, multiply terms that give rise to magnetic fields associated with 
dipole moments along Ha and perpendicular to it, respectively. In the limit e = Eq. ( p.4| ) 
simplifies [ |I| and the exact solution for $ is 

$ = + AiiifQliiOP^irj) cosif. (3.6) 

The parameter Ai\\ = A^ (for brevity we omit the upper index m = 1) is determined from 
the boundary condition = at ^ = ^o, the surface of the spheroid: 

^ ^ l + l/(l + eo)-eoarctan(l/eo)' ^^'^^ 

is always proportional to my, the longitudinal magnetic moment. For an ellipsoid, we 
show in Appendix A that 

m|| = ^fAi||. (3.8) 

It is instructive to verify that in the linear case the solution on Q leads to vanishing m±. 
From the linear relation between j and Vg in Eq. (|2.17| ) and the anisotropy in A as discussed 
II Q , the azimuthal dependence of Vg and j is identical. We can find the (/^-dependence 



m 



of the magnetic moment from Eq. ( |2.8|) and the appropriate element of integration given 
by Eq. ( |B6|) . The cp variable is separable. Consequently, the Jq^ dip integration can be 
performed analytically and yields m± = 0. Thus, any transverse component arises from the 
nonlinear terms. Their origin can be seen as follows: the superfluid velocity field has, in the 
presence of nonlinearities, an azimuthal dependence different from that in the linear case. 
This implies that the ip variable is no longer separable. The nonlinear terms in ( [A^ ) lead to 
higher harmonics for the (p dependence of Vg on Q. Then, it follows from Eq. ( p.4| ) that there 
is a small, but nonvanishing transverse field, H±, with a different azimuthal dependence from 
that in the linear case, which can contribute to m^. Therefore the nonlinear response of a 
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superconductor is responsible for Hx_ and, as a result, for m_L. In \ ri, the part of ifj, 
from which arises can be described by a transverse dipole, i.e. a potential of the form 
of the last term in Eq. ( |3.6| ), rotated by 90°: 

= A±/gi(^e)i^i(^7)sin</. (3.9) 

where in full analogy to the longitudinal case in Eq. (|3.8|) 

mx = ^fAix. (3.10) 

The potential $_l is the only contribution to from the general expression given by Eq. 
( p.4|) . This term is small, Ai_[_ <^ Ai\\ (recall the discussion about mj^/m\\ from subsection 
[11 B|) and does not exist when the penetration depths vanish, since the nonlinear effects are 
absent unless the field can penetrate the sample. Higher order multipole terms (/ > 1), as in 
( p.4|) do not contribute to the magnetic moment as can be seen from symmetry considerations 
or by explicit calculation using the orthogonality of P]^{ri). 

It is useful to recall explicitly the connection of the with the coefficients in a standard 
multipole expansion. In the spherical limit given by Eq. ( |3.3| ) these coefficients represent 
ordinary spherical multipoles. For large enough distances from Vt the asymptotic form of 
the / = 1 term is always a pure dipole. The spherical limit of Eq. ( p.6|) is 

$ = —HaT sin 9 cos ip H ^ smO cosip = —HaX H 5—, (3-11) 

with 

mo = -^Haa^ (3.12) 

where a is the sphere radius and we see once more that Ai\\{e = 0) is proportional to rriQ. 

For finite C, (not necessarily in the spherical limit), the magnetic moment of a spheroid is 
again obtained only from terms with / = 1 in the Eq. ( |3.4| ). The I = 1 terms, however, have 
spheroidal symmetry and their radial and angular dependence is not identical to -ij sin 6 as 
for a pure dipole in Eq. ( |3.11|) , they have an admixture of higher spherical multipoles. 
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At finite A the value for (proportional to my from ( |3.8| )) will be slightly different 
from Ai\\{e = 0), reflecting the difference between my and mo given by Eq. (|2.9a|) . The 
potential will also acquire additional terms as seen in Eq. ( p.4|) . The longitudinal terms 
with I > 2 and all of the transverse terms vanish at A = 0, that is, they are of higher order 
in e. They also vanish in the spherical case but only if A is isotropic and the relation j(vs) 
is purely linear. If, in addition, one includes the nonlinear terms in j(vs), so that the full 
relation (|3.4| ) applies, then the transverse dipole term appears even in a spherical geometry. 



IV. COMPUTATION OF THE MAGNETIC MOMENT 
A. General Considerations on Iteration Procedure 

The considerations from the previous section point to a method for obtaining a complete 
solution for m. This method, although it should work in principle, would naively lead to 
the need for an iteration which is in practice too cumbersome, and which we discuss here as 
motivation for introducing the perturbation method obviating the need for it. The important 
question here is that of calculating the magnetic dipole moment. In other boundary value 
problems to which our method can be extended, one must similarly determine the lowest 
nonvanishing multipole moments. 

To obtain a solution to Eq. ( p.6|) in the general case when the relation between j and Vg is 
nonlinear, given by Eq. ( 2.16|) (see also (|A2|) ) , one can in principle use the following iteration 



method. As a first step, one can solve these equations in Q with boundary conditions on 
dQ corresponding to the limit e = 0. That can be done by setting Ai\\ = Ai\\{e = 0) (and 
all the other j_ = 0), requiring continuity of the components if^ and at ,^ = ,^0) 
and also enforcing the condition on dQ, j ■ n = 0. These boundary conditions for the 
magnetic field are not exactly the desired ones, since continuity of the C, component cannot 
be demanded because the external field has been specified so that its ^ component vanishes 
at the boundary. The magnetic field outside the sample is simply obtained from Eq. (p. 2a ) 
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with $ given by ( |3.6| ). Inside one employs Eq. (|2.4]) to get the magnetic field from Vg. 
From the numerical solution in Q, obtained using the overrelaxation method discussed in 



IQ],^ we can compute the magnetic moment by using (p.8| ), and the fields from the above 
approximate solution. The computed value, m'^-'^^ the superscript (1) indicating the order 
of iteration, will be in general different from mQx: the computed magnetic moment is not 
the input value itlq. Hence, the actual problem has not been solved: the solution is not self- 
consistent. One can then imagine obtaining the correct solution from the following iteration 
process: Denote by {Aii\)^^\ (^ix)*-^-* the values of these quantities obtained from Eqs. 
( p.8|) , Eq. ( p.lO|) and the appropriate components of m^-'-^ Then use (Aiii)'^^^ (^ix)*-^-* in the 
computation of the exterior field, and use again this exterior field to solve Eq. ( p.6| ), repeating 
the procedure described in the previous paragraph. The second iteration yields e.g. 
(^i±)''^'' which in general would differ from (Aiy)'^^^ and {Ai_i_)^^\ Repeated iterations would 
give sequences {Aii\Y^\ and {Ai±Y^\... converging 

to the desired values of Ai\\, v4i_|_, when the moment generated by the computed currents 
equals the input value. At that point, the magnetic moment will be known. The higher 
order A's will not necessarily be known, but they do not contribute to m. In practice such 
a procedure could be implemented first for the larger longitudinal component and then for 
the smaller transverse component. 

B. Perturbation Method To Compute m 



The procedure described in [IV A| may be a lengthy and expensive process and it is for 
this reason that we develop, in this subsection, a procedure to bypass it. We compute my 
(accurate to first order in e) in a single iteration step, that is, a single pass through solving 
the equations inside the body as described above. 

The surface integral ( p.l3[ ), rather than ( |2.8| ) is the expression for m that is convenient 
for our purpose since as seen in subsection [11 B| , it divides m into two terms, mi and 1112, 
of different order in e. 
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It follows from the considerations of Section |I|, and it is the basis of our perturbation 
method, that to obtain m correctly to O(emo) it is sufficient to compute m2 from fields 
(i.e. H or Vg) accurate only to zeroth order. This is, as explained there, because a factor 
of e explicitly scales out of the expression for m2. Now, since the internal fields obtained 
by solving Eq. ( ^.61 ) at the first iteration level (as described in the previous subsection, i.e., 
with e = boundary conditions) are already accurate to the zeroth order, one iteration is 
sufficient to evaluate m2 at desired accuracy. The problem reduces, therefore, to that of 
correctly including the contribution mi to first order in e. 

To illustrate how this is done, let us recall the general form of the analytic solution for 
the spheroid. As seen in Section |T| (Eq. ( p.4|) ) it has the form of a multipole expansion 
with undetermined coefficients. Since the exact field H is continuous on dQ we can insert 



this general form in the expression ( |2.13|) for mi. Only the / = 1 terms, by virtue of ( p.8| 



and (p.lO|) , contribute to m, so that mi can be evaluated in terms of the unknown m. Thus 



we have 

m = mi(m) + m2 + O(e^mo), (4-1) 

where we emphasize that mi depends on the unknown / = 1 parameters, and where we 
introduce the overbar notation to denote quantities evaluated from the zero penetration 
limit external fields. Since terms of O(e^mo) can be neglected, it is possible to interchange 
m2 and m2 in the various expressions and we have done so. Expression ( [4. 1|) is an equation 
for the unknown m. 

To solve it in practice, consider the general form of $, Eq. ( p.4|) which indicates that $ 
(or equivalently H) can be separated into two parts: $ = + due to the applied field 
and to the presence of the superconductor, respectively. We can then write the contribution 
of these parts as 

mi(m) = mi($J + mi(<l>,). (4.2) 

We now define p by mi($a) = pmo. Since mo and mi($a) are now known, one can 
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determine the constant p which depends on the shape of fl, i.e. on the eccentricity. In the 
hmit e = 0, m2 = and from Eq. ([4.2|) we have the identity 



mo = mi = pmo + mo(l - p). (4.3) 

For e 7^ 0, when the solution for $ and H is given in terms of multipole expansion with 
unknown coefficients, mi($a) remains the same. The only difference in computing mi is 
that the coefficients in the terms arising from $r are now proportional to the correct, but 
still to be determined, value of the magnetic moment m, slightly changed from the e = 



case. All the remaining higher multipoles (/ > 1) of $, as mentioned in |T|, do not contribute 
to m and we have 

mi($r) = m(l -p). (4.4) 

Adding this to mi($a) we get 

mi = m — p(m — mo). (4-5) 
We can express Eq. ( ^.14| ) using ( ^.5) ) as 

m = m — p{m — mo) + m2, (4-6) 

so that we have the solution for m correct to 0(e), 

m = mo + -m2 ^ mo + -m2, (4.7) 
p p 

which determines all components of m. 

We illustrate this method using the textbook example of the isotropic, linear supercon- 
ducting sphere in a uniform applied field Ha, along the x-axis. In this case all the fields in 
Q satisfy the vector Helmholtz equation: 

V^F = ^F, (4.8) 

where F can be H, j, Vg. On the entire \ Q region, $r has a pure dipole form and H is 
given by taking the gradient of Eq. (|3.11| ): 
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2777- 

Hr = {Ha H ^) sin 9 cos (p, (4.9a) 



i/g = (Ha - ^) COS 0cos(^, (4.9b) 

TTl 

= {-Ha + ■^)sin<p, (4.9c) 

where r, ^, are spherical coordinates, and m a parameter to be determined. This is the 
general solution for any e (not just 0), in the field outside; there is only a dipole term in 
addition to that due to applied field. However, even if there were higher spherical multipoles 
in the general solution, that would not affect the evaluation of mi, since their contribution to 
m would vanish identically by symmetry. In the limit where the current does not penetrate 
into the superconductor, the magnetic moment is given by itiq = —\o?Ha (recall ( |3.12| )). 
Performing the elementary integral for mi in Eq. ( p.l4| ) we obtain 

12 1 
mi = -mo + -m = m- -(m-mo), (4.10) 

Comparing with ( [4. 51 ) we identify p = | in this case and using Eq. ( [4.7D 

m = mo + 3m2. (4-11) 

To determine the unknown m to 0(e) it remains to compute m2 and substitute it in ( [4.11D . 
Using the evaluation for m2, with the boundary conditions taken in the e = 0, from Appendix 
C we get the perturbation result for m: 

m = mo(l-3e), (4.12) 



where e = A/a. If we compare this to the expression for my given by Eq. ( |2.9a| ) we can 



read off a\\ = 3. This is the correct value for a sphere to this order, as given in textbooks [ 
P,P^. Thus, using the perturbation method with approximate boundary conditions in the 



evaluation of m2, we have obtained the correct value for m to first order in e. It is also 
instructive to calculate m analytically, with internal fields evaluated from ( p.6|) and e = 
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boundary conditions, Eq. ( p.l3| ). The calculation would yield m correct only to 0(mo); the 
value of ay is not correct; one gets ay = 2 instead of a;|| = 3. 

We return now to the oblate spheroid. In Q we allow the full nonlinear relation between j 
and Vg, given by Eq. ( [A2| ) . The magnetic field in R^\f2 that contributes to the computation 
of m is given in Appendix |B 3| . We aim to obtain the appropriate perturbation equation 
( [1.7D relating the unknown magnetic moment to m2, the term computed from the numerical 
solution in Q (using the first step in the iteration procedure, described in the previous 
subsection). To calculate the magnetic moment, we proceed as with the sphere example. 
We first evaluate the term mi from Eq. ( p.l4|) , recalling Eqs. (|3.7|) , ( |3l8|) , ( |3.10|) and n = ^. 



The integral, evaluated in spheroidal coordinates using (p6|), is elementary and we give only 
the result. The corresponding perturbation equation for m is, from ([4.7| ): 

m = mo + — ^m2, (4.13) 

or, writing its components explicitly: 

mil = mo + —^m2\\, (4.14a) 

m_|_ = m2^, (4.14b) 



P(6 



0) 

with 

p(0 = i(2 + - (1 + e)^arctanil/0), (4.15) 

p(^) is evaluated at the surface of the ellipsoid (^ = ^o)- As shown in Appendix |B 1| , is 
related to the eccentricity of the spheroid. In the spherical limit when ^ ^ oo we recover 
the spherical result, p = 1/3, obtained earlier. Another interesting limit is that of a flat disk 
(^ 0) where p{^) = 1/2. As discussed previously, mj_ — > for A — > 0. 

The above method applies not only to oblate spheroids, but to all geometries for which 
a general solution for $ in i?^ \ as an expansion in terms of orthogonal functions, (only 
one of them being dipolar at large distances) can be written. Furthermore, the method can 
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be extended, with one additional assumption, to a situation where the shape of Q precludes 
an analytic solution for the outside fields even at e = 0, and only a numerical solution, $, in 
that limit is available. We assume that, as in the analytic cases, the coordinate dependence 
of the terms in $ which contribute to the dipole m remains the same, up to a multiplicative 
constant, for e = and e 7^ 0. This assumption requires that the shape of Q produces no 
singularities in the fields. This might be a sufficient condition, but we know of no rigorous 
proof. 

The magnetic field at large distances, r ^ d, has the form given in Eq. tuq is the 



magnetic moment for Q, and all the higher multipoles can be neglected. We can again, as 
above (|4.2|) separate $ = + $r, where is the apphed field contribution. At r ci, 
$r is of dipolar form and the value of mo can be in principle numerically extracted either 
by using the left part of ( [4.3|) or from the asymptotic form. For e 7^ the potential $ has 
asymptotically the same dipolar form as $, but the unknown dipolar coefficient m differs 
from mQX. One can then proceed as in the analytic case. Consider as an illustration, the 
computation of my to 0(e). We can implement the method by writing the potential at e 7^ 
in the form: 

$(e ^ 0) = + ^(<l>,) + (4.16) 
mo 

where is a possible contribution to higher order multipoles only. Eq. ( [4. 161) merely 
expresses our assumptions in mathematical form. It can be better understood by recalling 
the discussion in Section |ITT| (and above in this Section) in particular the difference between 
Ai||(e = 0) and Aiy. 

Since mo is known, one can use $a to evaluate mi($a) (see (^^)) and hence the quantity 
p through mi($a) = prriQ. From the distribution of Vg on dQ we compute 777-2 and Eq. ( [4.7| ) 
gives the desired my. 

We see therefore that our method has considerable generality and our results can be 
summarized in terms of the following Theorem, the validity of which follows from the 
analysis in this section and the decomposition of m in [II B| . 
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Let us assume that: 

a) There exist a small parameter e ^ 1 and we consider Eq. ( |2.6| ) in Q that allows a 
sufficiently accurate solution in \ Q for ( |2.2b| ) with trivial Neumann boundary conditions 
on dQ, and at infinity — V$ = Ha, which satisfies the assumption discussed in connection 
with QOBI ). 

b) H in the interior of Q decreases with the distance from dQ not slower than exponential 
dependence given by a characteristic length ^ typical size of Q 

Then the following statements are true: 

l.It is possible to write m = mi + m2 as given by ( p.l4|) where mi is 0(mo) and m2 is 
O(emo). 

2. To obtain m from Eq. ( [4.7|) accurate to O(emo) it is sufficient to calculate the leading 
contribution to the term m2, the error in determining m being of O(e^mo). 

This perturbation method can be applied outside the field of superconductivity. For 
example, it is well known [ [I|] that an ordinary conductor in a high frequency, harmonic 
applied magnetic field (the frequency should satisfy quasi-static condition uj <^ c/d) behaves 
like a superconductor in a constant field. It is then possible to identify the small parameter, 
e <^ 1 as the ratio of skin depth, 6s, the typical length scale for field penetration in the 
conductor and the characteristic geometrical dimension, d. The computation of the magnetic 
field distribution is then achieved by solving the corresponding steady-state problem for a 
superconductor of the same shape, and m can be obtained using Eq. ( [4.7|) . 



V. NUMERICAL CONSIDERATIONS 

A. Dimensionless form of Equations on Q, 

In performing the calculations and describing the results, it is convenient to introduce 
dimensionless quantities. We recall that f2 is a fiat spheroid and with the magnetic field 
applied in the x — y {a — b) plane, most of the current will fiow parallel to the a — b plane 
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and its decay will be governed by Xab- It is therefore convenient to measure the length in 
the units of Xab- We then define dimensionless fields V, J, and H: 

Vc Ho AnXab 

where Vc = Ao/f/ is the critical velocity (discussed in Appendix A), Aq is the amplitude 
of the order parameter, defined in Appendix A, and we have introduced a characteristic 
magnetic field Hq as 

H„ = (5.2) 

Xabkab 

where 0o = ttc/i/c is the flux quantum and S,ab = Vf/nAo is the in-plane superconducting 
coherence length (not to be confused with ^, the spheroidal coordinate). The definition ( p.2| ) 
involves precisely the same numerical factors as that used in [ Q . The required equations 
are easily rewritten in terms of these quantities. Equation ( p.6| ), using the relation between 
j and Vg given by ( |A^) , then becomes 



(V X V X V)^./y = -Kr/y(l - ti = -Vx' ,y' + iV^'y, (5.3a) 

((5V X V X V), = -v;(i - i ) = + 

I Vx' \ + \ Vy'\ 

where 5 = {X(./XabY — ^c/^ab and we define Nxi-, Ny/, as the terms nonlinear in the 
velocity. The equations are written in terms of the primed coordinates and the derivatives 
are with respect to the dimensionless length measured in units of Xab- 



Before discretizing Eqs. ( |5.3|) we transform them to oblate spheroidal coordinates. We 
start by writing these equations in the unprimed {x, y, z) coordinate system where, we recall, 
the X-axis lies along the applied field. The linear part of the equations looks identical in 
primed or unprimed coordinates and we only need to carefully transform Nx\y',z to A^^;, Ny 
and Nz terms nonlinear in the velocity along the unit vectors x, y and z respectively. We 
have 

Nx = ti{Vx' IK'I cos* + Vy' |V^y I sin*), (5.4a) 
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Ny = IK'I sin^ - Vy> IK^'I cos^), 



(5.4b) 



= t2V,- 



V3 + Vy'^ 



\v,'\ + \Vy,y 

where = VxCos"^ — V^sm\l/, Vy' = V^sin\l/ + V^cos^l/ and V^i 
the components of velocity in spheroidal coordinates: 



(5.4c) 



Vz- Or, if we express 



= cos(v3 + ^)(a\4 + dVr,) - sin(v? + ^)K^ 



(5.5a) 



Vy> = sm{ip + ^)(a14 + dVr,) + cos(y? + ^)V"^, 



(5.5b) 



Vz = -dV. + aVr 



(5.5c) 



where a = ^^^^aff''" and ci ^ 



We can now write the nonlinear part (from equation 



]6D and (|5.4| )) along each spheroidal coordinate. For example, along ^ we get 



= a(cos ifNrc + sin fNy) — dN^. 



(5.6) 



^x,y,z are entirely expressed in terms of spheroidal components as shown above. In an 
analogous way we can obtain the remaining nonlinear parts Nr^,,p- 

Using Eqs. ( |3.2| ) we transform the inverse of the penetration depth tensor given in 
Cartesian coordinates by a diagonal tensor with components (A^ff , A^^^, A~^), (recall Eq. 
( p. 171 )) to spheroidal coordinates: 



A 



-1 



Pi P2 
P2 P3 
1 



where pi,2,3 are defined by 



Pi 



P2 



^2 + ^2 



/2 



(5.8a) 



(5.8b) 
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PS = " ■ (5-8C) 



The resulting form of the dimensionless equations in spheroidal coordinates (which we will 
solve numerically) is 

/2(V X V X V)5 = -fiipiV^ + P2V,) - N^), (5.9a) 
/2(V X V X V), = -f{{p2n + PsK,) - N,), (5-9b) 

/2(V X V X V)^ = -f\V^ - N^). (5.9c) 

Expressions for the differential operator /^Vx Vx V in spheroidal coordinates are included in 
Appendix P 2| . Equations ( |5.9|) have to supplemented with appropriate boundary conditions, 
as discussed in Section II. The boundary condition at infinity is satisfied by the use of the 
analytic solution for R'^ \ Q. Since Xab <^ C we can put V = at ^ = 0, as all the fields 
vanish deep inside the sample. Continuity of the Ti field on dQ is achieved through 

vxv = n\^=^,, (5.10) 

where the right hand side is the external dimensionless field at the surface of the ellipsoid. 
The remaining boundary condition J • n = = on dQ is generally nonlinear and readily 
enforced by observing that the RHS of Eq. ( |5.9aD is oc J^. 



B. Computational Grid and Discrete Variables 

The implementation of the perturbation method from Section ^ is not restricted to a 
particular algorithm for solving the relevant equations in Q. In the remaining part of this 
section we outline as an example one suitable algorithm using a modification of the Gauss- 
Seidel relaxation method. We first discuss the discretization of the complete nonlinear, 
three-dimensional (3D) problem. It is then possible to consider, as a special case, the two- 
variable discretization of the linear problem were the ip dependence is known analytically. 
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The numerical solution to such a problem is then used as the initial guess for the relaxation 
method of the complete 3D problem. 

Eqs. ( |5.9| ) or their counterparts ( |2.6| ) and have definite parity: Vx-, Vy are even and 

Vz is odd in z. It is therefore sufficient to consider only the upper, positive z or the 

lower, negative z half of Q and extend by parity the obtained solution to the whole 
Q. The computational domain G is obtained by parameterizing the physical domain f2_ 
using oblate spheroidal coordinates (^, t], ip). We consider a uniform grid on G with mesh 
widths /ig, hj^ and h^p. The choice of a grid uniformly spaced in variable 77, generates denser 
grid points corresponding to the part of Vt^ (the r/ ^ region) with higher curvature and 
greater field variation. An arbitrary grid point on G is given by (^j, ?7j, ipk) or just (i, j, k) 
for brevity: 

y^ij,k = iii + 'n3fl + Vk0, (5.11) 

where the grid coordinates are given by 

ii = ih^, T]j = -1 + {^+j)hr^, (pk = --^ + kh^, (5.12) 

and the indices run through values i = 0, ..,n^, j = 0, ..,n^ and k = 0, ..,n^. Mesh widths 
are given by = ^, /i„ = ^ and = J. 

The discrete variables are denoted by the same symbol as their continuous counterparts, 
for example, V^.j^ fc represents at the grid point (^j, rjj, ip^). The discretized approxi- 
mations of derivatives used have second order accuracy in the mesh widths. We will use 
the letter D to represent the discretized approximation, upper indices 0, +, — denote the 
central, forward and, backward approximation respectively, and the lower indices denote 
the corresponding variables of differentiation. For example, denotes the mixed differ- 
entiation with respect to ^ and 17, where the central difference formula is used in ^ and the 
forward difference formula in rj variable. For the interior grid points we only use the central 
difference formula for all the derivatives and omit the upper indices. We also use central 
differences for derivatives at all the grid points on dG that do not require introduction of 
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fictitious grid points ^ G. For example, in computing at (n^, j, k) we would need to use 
a grid point at z = ra^ + 1 ^ G, we avoid that by using a backward difference at {n^,j, k). 
Similarly we employ forward differences were necessary. 



C. Implementation of Boundary Conditions and Equations on G 

The grid boundary, dG, consists of six two-dimensional planar surfaces with grid points 
described by (0,j, /c), {n^J,k), (i,0, /c), {i,nr„k), {ij,0) and (i,j,n<^). 

On the surface ^ = 0, {0,j,k), which corresponds to the equatorial {z = 0) disk in fl 
with radius equal to the focal length /, we impose trivial Dirichlet boundary conditions. 
As we have discussed in Section II, this follows from the requirement that deep inside Q all 
fields should vanish. This eliminates possible difficulties from the singularities of the various 
differential operators at = 0. Any remaining singularities of the equations on G would 
come from points with coordinates rj = and rj = ±1. On the surface rj = {i, j = rirj, k), 
corresponding to part of the z = plane, we have from the known parity of V: 

V^^^ = 0|^=o, (5.13a) 

d^V^ = 0|^=o. (5.13b) 
We implement Eq. (^.13b| ) as D:^Vr,-ij=n^^k = 0: 

and in the iterative solution we write down explicitly K?;jj,A; from this expression. The 
remaining region on G that could result in singularities of the differential operators is at 
T) = —1, which corresponds to a line through the "south pole" and the origin of Q, i.e., a 
segment starting at the origin and ending at a point with Cartesian coordinates (0,0, z). 
The choice of grid given in Eq. ( |5.12| ) excludes this segment, the closest point on the grid 



is hfj/2 away. The numerical solution for V in the vicinity of 77 = —1 is well behaved (as 
it is for the linearized equations in the geometries that permit an analytic solution in Q). 
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It is therefore possible to extrapolate the obtained numerical solution to rj = —1. At the 
grid boundary surface = 0, k) we use the forward difference formula for the derivatives 
with respect to t], is obtained analogous to Eq. (|5.14| ) by replacing /i^ with — /i^ and the 
indices j — 1, j — 2 by j + 1 and j + 2 respectively. The second derivative D^^ is taken as: 

where F represents any component of a vector field. We can obtain similar formulae for D^^ 
and D^^. At the surface boundary k = 0), where (p = — vr, we proceed in an analogous 
way, the derivatives with respect to (f are expressed with forward differences. The other 
part of dG with ip = const i.e., {i,j,k = n^p) corresponds to the same boundary surface 
[ip = —tt) and we can impose the simple periodic boundary conditions 

= i^„-,o. (5.16) 
On the remaining part of dG, the surface {i = n^,j, k) aX = ^o? we impose continuity (as 



discussed in section IV) of the rj and (p components of 7i. From Eq. ( [5.10| ) we can express 
Vrj-ij^k and V^p-ij^k respectively to obtain their updated values in each step of the relaxation 
procedure. The equation for V^-ij^k is obtained from J • n=0 by setting the RHS of Eq. ( [5.9aj ) 
to zero 

PiV^ + P2Vrj - = 0. (5.17) 



In the relaxation procedure, the nonlinear term A^^ is included using Picard's method [ ^ 
the value of is taken from the previous iteration. If we denote by an upper index n the 
number of the iteration (in the relaxation procedure), the nonlinear boundary condition can 
be implemented as 

Pi Pi 

The nonlinear and nonanalytic terms can be simply included by using Picard's method. In 
addition to the nonlinear boundary condition ( |5.18| ) we shall also use Picard's method to 
include the nonlinearities stemming from Eq. (|5.9|) . 
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For the grid points ^ dG it is possible to use central differences. In the LHS of Eq. (^.91), 
given explicitly in Appendix C, we replace each partial derivative by the appropriate central 
difference, D. 



D. Modified Gauss-Seidel Relaxation 

The solution to the nonlinear, 3D problem using the relaxation method consists of two 
steps. The first is obtaining a good initial guess using the numerical solution to the linear 
equations given by Eq. (|2.6| ) and ( p . 1 7| ) . The ip dependence is then known and it can be 



separated out. The method of successive overrelaxation can be applied to the resulting two- 
variable problem in the coordinates ^ and rj. The boundary conditions on the continuity of 
Ti, as discussed in the previous subsection, are implemented: the magnetic field outside is 
taken in the e = limit {Ai\\ = y4i||(0)). In the relaxation procedure each component of V 
is expressed in terms of the corresponding component of the linearized {N^^^^^p = 0), two- 
dimensional form of Eq. (p.9D so that the resulting matrix equation is diagonally dominant 



2^ . For a detailed discussion of the method see [ The numerical solution, the 

distribution of V{^i, rjj), is supplemented with the known ip dependence and my is calculated 
using the perturbation method. 

The second part of the algorithm is a modification of the Gauss-Seidel (the overrelaxation 
parameter, uj = 1) method for the full 3D problem. The previously obtained solution for the 
linearized equations is the initial guess for Vij^k, and the value m\\ i.e., the corresponding Ai\\, 
is used in the expression for the Ti outside. If we denote the linear part of Eq. (|5.9|) as L(V) 
and the nonlinear part as N(V) the relaxation procedure can be described symbolically as: 

L(V("+')) = N(V("^), n=l,2,.. (5.19) 

After completion of each relaxation step we update the old values, V*^"\ at each grid point as 
N(V^"'') = N(V*'"^^^). As done with the boundary condition from Eq. ( [5.18| ) we use Picard's 



method to include the nonlinear terms which are also nondifferentiable. The numerical 
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solution for the linearized equations is a very good initial guess for the full problem, since 
the nonlinear terms are small, which compensates for the relatively slow convergence of 
Picard's method. 

VI. NUMERICAL TESTS 

For our computations we used the Cray C90 of the Minnesota Supercomputer Institute. 
We first tested the computer code on the example of the linear, isotropic sphere with applied 
field along the x axis, where the analytic solution is known. The spherical geometry was 
realized as the spherical limit of the spheroid. We used = 1000 and / = 0.1, so that the 
radius of sphere was a — > /,^o = 100 (in units of A), the corresponding e = X/a = 10~^ 
and the ratio of the spheroidal semi-axes is A/C = 100.005/100. We used the two-variable 
version of the code in the variables ^ and rj. The computational grid spanned a spherical 
shell of thickness 7 (A), and because fields decay exponentially away from the surface, we 
imposed trivial Dirichlet conditions at ,^ = 930. We used = 200 and = 50 for the 
number of grid points. To test the convergence of the relaxation algorithm, we tried the 
very poor initial guess of zero fields everywhere on G. The boundary conditions were taken 
in the e = limit, that is Ai\\ = Ai\\{e = 0) from Eq. ( p.7|) . The overrelaxation parameter 
was uj = 1.8 and after 1000 relaxation steps (43 sec of CPU time) we obtained m using Eq. 
( p.l4| ) and the perturbation method with Eq. ( [4.11|) . In the term m2 the distribution of 
Vrj, K/p at ^ = ^0 (from the numerical solution on G) was extended by parity to the entire dfl 
and supplemented with the known ip dependence. Thus, integration over (p was performed 
analytically and that with respect to rj, numerically. The extracted constant (from Eq. 
( [4. lip ) was a\\ = 3.1, within 3% of the analytically obtained value of 3 from Eq. ( |4.12| ). The 
use of a conventional procedure to compute m would require repeating the whole iteration 
procedure to get an improved value for a\\, as we have in principle described in [IV B . 

We also verified that the numerical solution for the full three-variable problem in this 
geometry and exact boundary conditions, has the correct form. The angular dependence 
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of the solution was J^, oc rj simp and J^,, oc (1 — rj'^Y^'^ cos ip which in the spherical 
limit T] —>■ cos 9 corresponds to the analytical solution for a sphere. The numerical solution, 
for the range of grids considered, was accurate between three and four significant figures for 
every grid point where J (V) was numerically significant. 

For the physical results related to the nonlinear response of an oblate spheroid we refer 
to [ ^, and discuss here only some aspects not covered there, that illustrate the numerics. 
Analysis of the transverse magnetic moment shows that [ ^ 

cx ^H,F{m), (6.1) 

from which we infer that a±_ oc Ha/Ho where is the angular dependence on \& 

which has [ Q period 7r/2. We will give results for \E' = vr/S, approximately the maximum 
of F. The longitudinal moment, due to field penetration, differs from tuq and it can be 
characterized, as we have shown earlier, by the parameter a\\. This parameter includes 
contributions from the linear part of j(v), independent of Ha, and from the nonlinear part, 
dependent both on Ha/Ho and on \E'. 

We consider an oblate spheroid with = 0.144338 {A/C = 7) at Ha/Ho = 0.1 (in the 
experimentally relevant range) and \E' = 7r/8. For / = 1000 (in units of Xat, defined in [11 (J 
and |V A|) we have used = 550, = 50 and = 30. The results for a\\ and a± are 
given as a function of the material parameter 6 = {Xc/^abY in Table I. For fixed .^o (fixed 
shape), we have considered various sizes of spheroid (different /), changing e = Xab/C ^ 1 
by a factor of four (at fixed Aaf,, Ac) and verified that a\\, a± are size independent within 
numerical accuracy. 

In the next two figures we display some of the numerical results for the field distributions 
calculated under the same conditions and with the same parameter values as in the previous 
paragraph. In Fig. 3, we show results for the current at surface of the spheroid, (^ = .^o)- 
The components J^,r]{^o, Vi at the fixed azimuthal angle ip = 48°, are shown as functions 
of rj. For comparison we recall here also the corresponding angular behavior for a sphere, 
as given earlier in this Section. In Fig. 4 we show some results for the magnetic field as it 
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penetrates into the sample. We plot the two components TCrj^^iC, rj, ip) as functions of — ^ 
at constant ip = 0°, r] = —0.693, (see Fig. 2). The plot illustrates the difference between 
the components of the field arising purely from the linear equations and those which are 
due to the nonlinear effects. The component (at ip = 0°) is very predominantly "linear" , 
and displays exponential-like behavior. The other component plotted, vanishes in the 
linear case (at ip = 0°), and it arises solely from the nonlinear effects. This behavior is far 
from being an exponential; its derivative changes sign, for the same physical reasons as the 
nonlinear current does, as discussed in [ § . 

VII. CONCLUSIONS AND FUTURE WORK 

In this paper, as the main result, we have presented a perturbation method to compute 
the magnetic moment of a bulk nonlinear and anisotropic superconductor. This method 
could be implemented in conjunction with various algorithms for solving boundary value 
problems in electrodynamics. Suitable generalizations would certainly increase the range of 
its applicability from that discussed in this paper. Obvious examples include considering in 
detail other shapes of Q and computing higher order multipoles. We have showed that our 
method increases the accuracy of computation while very significantly reducing the required 
computational work. 

In this paper the numerical example of an oblate spheroidal geometry was discussed in 
detail to illustrate the perturbation method and also to give guidelines for possible improve- 
ments. The numerical algorithm which was employed for solving the nonlinear Maxwell- 
London equations provided more than sufficient accuracy, as seen from experimental con- 
siderations: the uncertainty of the input experimental parameters significantly exceeds the 
accuracy of the results obtained. For our computations we have used Cray C90 and memory 
requirements were not the limiting factor. It is however possible to make various improve- 
ments to the numerical algorithm. One could consider nonuniformly spaced grid points 
along the ^ coordinate, denser close to the ^ = ^o, in order to reduce the overall number 
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of grid points. It might also be advantageous to replace the Gauss-Seidel relaxation on G 



by some other method, such as GMRES [ pH] or one of the various multigrid algorithms 



|T9|. In future work we will consider some of these improvements and investigate possible 
generalizations of the perturbation method presented in this paper. 
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APPENDIX A: THE NONLINEAR RELATION j(vs) 



In order to express j(vs) from Eq. ( p.l6|) , we introduce two coordinate systems: x — y, 
fixed in space such that the applied field is along the x-axis, and x' — y' which is fixed to 
the crystal. We consider an order parameter of the so called (i-wave form: 

A = Aosin(20), (Al) 

where is the azimuthal angle referred to a node and Aq the gap amplitude. It has been 
shown [ ^ that in the field range of experimental interest and with suitable assumptions for 
the Fermi surface, Eq. (|2.16| ) can be rewritten, at sufficiently low temperatures as 



jx',y' = -epabVx',y'il - — \Vx',y'\), (A2a) 



V, 



c 



j^ = -epMl-- "\lf X (A2b) 

Vc \Vx' \ + \Vy'\ 

where ti and ^2 are constants, ti = |f , ^2 = ||- The critical velocity is Vc = Aq/vj. We 

2 _„ _\2 

define here pat = ^>^ab and pc = pab^- 
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APPENDIX B: QUANTITIES IN SPHEROIDAL COORDINATES 



1. Magnetic Moment 

We show here that the expression for the longitudinal magnetic moment for an oblate 
spheroid, given by Eq. ( p.8| ), is equivalent to the standard result expressed in terms of the 
demagnetization factors. For a field applied along the x axis, in the limit A = 0, we have [ 

moi, = mo. = -4 (Bl) 



''X J 



where V is the volume of a spheroid and 

1 1 + 6^ 

nx = ^^^^ ~ (^rctan{e)) (B2) 

is the appropriate demagnetization factor. The eccentricity is given by e = [A^/C"^ — 
= 1/^0- In terms of spheroidal coordinates we have A = /(I + ^qY^"^, C = f^o and 
V = ^P{1 + ^q)C,o- Including these expressions in Eq. (pl| ) we get 



3^'^" 1 + 1/(1 + ^0') -Coarctan{l/^o) ^^^^ 



and we recover my = ^f^Ai. 



2. Integral and Differential Operators 

The metric coefficients in oblate spheroidal coordinates are given by [ |T7|: 



9u = fYT^' ^^^""^ 



922 = r-r^. (B4b) 
1 — r/^ 

9z, = f{l+e){l-v')- (B4c) 
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One can then calculate the appropriate elements of integration and differential operators. 
For example, 



/ ds= f dr^ r d^ni + ef'\e + v'f'\ m 

Jan J-i Jo 



f dn = f'dU dr^r dipf^e + V^)- (B6) 
Jn Jo J-i Jo 

$,0 corresponds to the value of ^ at the boundary dQ. The gradient is 

^ - /(^2 + ^ + /(^2 + ^2)1/2 W', + ^ ^2)1/2(1 _ ^2)1/2 '^^'^^ ^^^^ 



To solve equations (2^) the expression V x V x v should be transformed to oblate 
spheroidal coordinates: 

/^(V X V X v)^ = aodr,r,v^ + aidrjV^ + a2d^^v^ + a^v^ + a^d^^r^Vr, + a^d^Vr, + aQdr,Vn 

+ a-jVr^ + asd^^v^ + agd^v^, (B8) 

/^(V X V X v)^ = bod^^v^ + + &2f^,7^^5 + ^3^5 + hd^^Vr, + hd^Vr, + hd^^Vr^ 

+ foyt;^ + bsdn^v^ + bgd^v^, (B9) 



/^(V X V X v),^ = Pod^^v^ + pid^v^ + P2dn^Vn + Psd^v^ + p^d^^v^ + p^d, 



+ Pedirv^ + prdr^v^ + ps-^.^- (BIO) 

We recall that / is the focal length scale factor. The coefficients Oj, bi,pi are given by (using 
the abbreviations u = {1 + ^^)^/^ s = (1 - r/^)^/^ w = {C,'^ + rff/"^): 



ao = -^^4 = tbo = -^64 = ^2^4 = P5 = -f;2 , 
ai = -jbr, = 7P6 = -P7 = -5, 



^2 = -liag = 



1 



2„2 ; 



_ 2-/7^-g2+3g2^2 

03 - ^;;6 , 
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5g = ivs(3+2e-V^) 



h 



w 



6 



Pi = — P3 = 1-3 



(Bll) 



3. Magnetic field which contributes to m 

We write down here explicitly the contributions to the magnetic field which arise from 



the potential $ given by the sum of ( p.6|) and (|3.9| ). From Eq. ( 2.2a ) and ( P7|) we have 



^€ = cos^ + h±{0 sin^), (B12a) 

Hv = - (^2/^2) 1/2 (/^(O cos^ + hMO sin^), (B12b) 

= - ^^2)1/2 (/^(O sin^ - cos^), (B12c) 
where the functions fi{C), /2(0 given by: 

= H,^ + + -^7^ - ^arctan{l/0), (B13a) 

= + - earcton(l/e)), (B13b) 

/2(0 = H^{1 + + (T^fiyiT^ - (1 + e^)^/'«rctan(l/0), (B13c) 

/2x(0 = - (1 + ef'^arctan{l/i)). (B13d) 

Terms with /i and /2 represent the longitudinal part of H and those with and /2± the 
transverse part. 
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APPENDIX C: EVALUATION OF ma TO O(emo) FOR A SPHERE 



We evaluate here the magnitude rri2 for the example of an isotropic spherical super- 
conductor with a linear j(vs) relation. We take the field along the z direction, since the 
magnetic moment m2 is independent of this choice. As explained in the text, we must find 
the internal fields by solving ( |2.6| ) (which in this case reduces to the Helmholtz equation) 
with boundary conditions enforcing continuity of and Hq and external fields calculated 
in the e = limit. By symmetry, if^ = and we need only to impose the continuity of Hg 
at the boundary, r = a. We obtain: 

\ 2 2 

He = -A—{{1 + sinh(r/A) - ^ cosh(r/A)) sm9 = {-Ha + ^) sin^. (CI) 



A is a constant to be determined and rriQ is given by Eq. (|3.12|) . Since e = ^ <^ 1 we can 



approximate sinh(a/A) ~ cosh(a/A) ~ |e''"/^\ Keeping only the leading term in the LHS of 



A = SiY^e-^'^/^) (C2) 



Eq. (|C|) we get 



and 



We can now compute 



Hg = -lHa-e-^^-^^/'sme. (C3) 
2 r 



m2 = — / dQH (C4) 



57r Jn 



Integration is performed in spherical coordinates giving 

m2 = ^Xa^Ha = ~~"^o ~ O(emo). (C5) 

This result (O(^mo)) was expected since for exponentially decaying fields, integration over 
the whole volume of Q is effectively only integration over the region ~ A away from its 
surface. 
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FIGURES 

FIG. 1. Geometry considered here. The superconducting region is an oblate elhpsoid of 
revolution. The x, y and z directions are fixed in space. The field is applied along the x axis, as 
indicated, while is along the y axis. The long and short semiaxes values are called A and C in 
the text, respectively. 

FIG. 2. The relation between oblate spheroidal rj, ip) and Cartesian coordinates is illustrated 
at a fixed azimuthal angle (p = 0°. The surfaces ^ = constant are oblate ellipsoids of revolution 
around the z-axis with semiaxes A = /(I + ^^)^/^, C = f(,. The surfaces |r/| = constant are 
one-sheeted hyperboloids of revolution. The surfaces ip = constant are planes through the z-axis 
making an angle ip with the x — z plane. 

FIG. 3. The spheroidal components of the dimensionless current (see ( ^.l]) ) at the surface of the 
spheroid = ^q)- The quantities plotted are J^(^0)^?)¥' = 48°) (solid line) and JipiCo^V^V = 48°) 
(dashed line). We have used H/Hq = 0.1, A/C = 7,^ = vr/S, 5 = 16. 

FIG. 4. Illustration of the nonlinear effects on the magnetic field. The solid line is the compo- 
nent firfiCv = — 0.693, = 0°), which is very predominantly "linear", normalized to its maximum 
value, while the dashed line is 'HipiCv = — 0.693, = 0°), which arises solely from nonlinear ef- 
fects, also normalized to its own, much smaller, maximum value. Both components are plotted 
as functions of D = — The non-exponential behavior of the dashed line is a signature of its 
nonlinear character. All parameters used are as in the previous Figure. 
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TABLES 

TABLE L The parameters ay, a±_ which determine the magnetic moment for e ^ (see Eq. 
(|2.9| ), given as functions of the material parameter 5 = {Xc/Xab)'^, computed for an oblate spheroid 



with A/C = 7 at a field H/Hq = 0.1. 

5 = {K/Kh? ay lO^a^ 

16 L9 3.8 

25 2.1 3.6 

36 2.2 3.3 

50 2.3 3.1 
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